iSSF Vegetation Model Plotting
Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain vegetation categories as a covariate, with interactions for movement parameters.
Load packages
Import data
New names:
Rows: 1062908 Columns: 42
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(4): pheno_start, pheno_end, season, day_period dbl (35): ...1, id, burst_,
x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr... lgl (1): case_ dttm (2):
t1_, t2_
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`Prepare the data
Code
# scale the covariates
ssf_dat <- ssf_dat %>% mutate(srtm_start_scaled = scale(srtm_start),
srtm_end_scaled = scale(srtm_end),
water_dist_start_scaled = scale(water_dist_start),
water_dist_end_scaled = scale(water_dist_end),
ndvi_start_scaled = scale(ndvi_start),
ndvi_end_scaled = scale(ndvi_end),
nbr_start_scaled = scale(nbr_start),
nbr_end_scaled = scale(nbr_end),
sl_scaled = scale(sl_),
log_sl_scaled = scale(log_sl_),
cos_ta_scaled = scale(cos_ta_)
)
# get scaling factors
sl_scale_sd <- attr(ssf_dat$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_dat$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_dat$cos_ta_scaled, "scaled:scale")
# rename the pheno categories to set a different reference level
ssf_dat$pheno_end <- relevel(as.factor(ssf_dat$pheno_end), ref = "shrub_grass")
ssf_dat$pheno_start <- relevel(as.factor(ssf_dat$pheno_start), ref = "shrub_grass")
unique(ssf_dat$pheno_start)[1] shrub_grass dry_grassland tree_grass bare_ground open_forest
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass[1] shrub_grass dry_grassland tree_grass bare_ground open_forest
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grassPlot the used-available distributions
Code
# checking for all individuals
buffalo_ids <- unique(ssf_dat$id)
# inspect the data for dry season
ssfdat_dry %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "Case")+
ggtitle("Dry season", subtitle = "All buffalo")+
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Inspect the data for each buffalo
Code
for(i in 1:length(buffalo_ids)) {
print(ssfdat_dry %>% filter(id == buffalo_ids[i]) %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "case_")+
ggtitle("Dry season", subtitle = paste0("Buffalo ", buffalo_ids[i])) +
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
)
}`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Code
ssfdat_wet %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "case_")+
ggtitle("Wet season", subtitle = "All buffalo")+
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
For each individual
Code
for(i in 1:length(buffalo_ids)) {
print(ssfdat_wet %>% filter(id == buffalo_ids[i]) %>%
group_by(case_, pheno_end) %>%
summarize(n = n()) %>%
mutate(prop = n / sum(n),
label = paste0(round(prop * 100, 1), "%")) %>%
ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) +
geom_col(position = position_dodge2()) +
geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
labs(x = "Land use class", y = "Proportion", fill = "case_")+
ggtitle("Wet season", subtitle = paste0("Buffalo ", buffalo_ids[i]))+
scale_fill_brewer(palette = "Paired", name="Case",
breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
theme_light() +
theme(axis.text.x = element_text(angle = 20, hjust = 1))
)
}`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.
Read in the saved model objects
Model summaries
Model fitting - interacting with movement
Code
tic()
buffalo.tmp <- glmmTMB(case_ ~
pheno_end +
sl_scaled +
log_sl_scaled +
cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id) +
# (0 + pheno_end | id),
diag(0 + pheno_end | id),
family = poisson(),
data = ssfdat_dry,
doFit=FALSE,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA,1:(nvarparm-1))))
# to render the script we comment out the model fitting, as we are loading a saved model
buffalo.pheno.dry.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)
summary(buffalo.pheno.dry.mvmt) Family: poisson ( log )
Formula:
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id) + diag(0 + pheno_end | id)
Data: ssfdat_dry
AIC BIC logLik -2*log(L) df.resid
1160974.1 1161257.6 -580462.1 1160924.1 620298
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
step_id (Intercept) 1.000e+06 1000.0000
id pheno_endshrub_grass 4.206e-02 0.2051
pheno_endbare_ground 1.289e-02 0.1136 0.00
pheno_enddry_grassland 2.809e-02 0.1676 0.00 0.00
pheno_endopen_forest 1.762e-01 0.4197 0.00 0.00 0.00
pheno_endtree_grass 2.464e-02 0.1570 0.00 0.00 0.00 0.00
Number of obs: 620323, groups: step_id, 56393; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.905077 4.211376 -0.690 0.490310
pheno_endbare_ground 0.289401 0.064373 4.496 6.93e-06
pheno_enddry_grassland 0.026845 0.071515 0.375 0.707379
pheno_endopen_forest -0.023060 0.123586 -0.187 0.851983
pheno_endtree_grass -0.011854 0.070324 -0.169 0.866140
sl_scaled 0.017248 0.007990 2.159 0.030861
log_sl_scaled -0.003373 0.008243 -0.409 0.682367
cos_ta_scaled -0.001107 0.006586 -0.168 0.866556
sl_scaled:pheno_startbare_ground -0.024199 0.017730 -1.365 0.172290
sl_scaled:pheno_startdry_grassland 0.058082 0.016770 3.464 0.000533
sl_scaled:pheno_startopen_forest -0.033550 0.016937 -1.981 0.047607
sl_scaled:pheno_starttree_grass 0.018553 0.018620 0.996 0.319039
log_sl_scaled:pheno_startbare_ground -0.040732 0.015809 -2.577 0.009981
log_sl_scaled:pheno_startdry_grassland -0.094216 0.017167 -5.488 4.06e-08
log_sl_scaled:pheno_startopen_forest 0.074716 0.017826 4.191 2.77e-05
log_sl_scaled:pheno_starttree_grass 0.013739 0.019710 0.697 0.485745
cos_ta_scaled:pheno_startbare_ground -0.046518 0.012831 -3.625 0.000289
cos_ta_scaled:pheno_startdry_grassland 0.009236 0.013844 0.667 0.504667
cos_ta_scaled:pheno_startopen_forest -0.038082 0.013482 -2.825 0.004732
cos_ta_scaled:pheno_starttree_grass -0.035782 0.014887 -2.404 0.016233
(Intercept)
pheno_endbare_ground ***
pheno_enddry_grassland
pheno_endopen_forest
pheno_endtree_grass
sl_scaled *
log_sl_scaled
cos_ta_scaled
sl_scaled:pheno_startbare_ground
sl_scaled:pheno_startdry_grassland ***
sl_scaled:pheno_startopen_forest *
sl_scaled:pheno_starttree_grass
log_sl_scaled:pheno_startbare_ground **
log_sl_scaled:pheno_startdry_grassland ***
log_sl_scaled:pheno_startopen_forest ***
log_sl_scaled:pheno_starttree_grass
cos_ta_scaled:pheno_startbare_ground ***
cos_ta_scaled:pheno_startdry_grassland
cos_ta_scaled:pheno_startopen_forest **
cos_ta_scaled:pheno_starttree_grass *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1130.739 sec elapsedCode
Code
[1] "Model location: buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_dry_2025-05-20.rds"Code
tic()
buffalo.tmp <- glmmTMB(case_ ~
pheno_end +
sl_scaled +
log_sl_scaled +
cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id) +
# (0 + pheno_end | id),
diag(0 + pheno_end | id),
family=poisson(),
data = ssfdat_wet,
doFit=FALSE,
control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA,1:(nvarparm-1))))
# to render the script we comment out the model fitting, as we are loading a saved model
buffalo.pheno.wet.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)
summary(buffalo.pheno.wet.mvmt) Family: poisson ( log )
Formula:
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +
pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
(1 | step_id) + diag(0 + pheno_end | id)
Data: ssfdat_wet
AIC BIC logLik -2*log(L) df.resid
827189.3 827464.3 -413569.6 827139.3 442560
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
step_id (Intercept) 1.000e+06 1.000e+03
id pheno_endshrub_grass 2.569e-02 1.603e-01
pheno_endbare_ground 7.612e-02 2.759e-01 0.00
pheno_enddry_grassland 2.976e-05 5.455e-03 0.00 0.00
pheno_endopen_forest 7.819e-02 2.796e-01 0.00 0.00 0.00
pheno_endtree_grass 1.721e-02 1.312e-01 0.00 0.00 0.00 0.00
Number of obs: 442585, groups: step_id, 40235; id, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.142132 4.985576 -0.630 0.528535
pheno_endbare_ground 0.500199 0.087277 5.731 9.97e-09
pheno_enddry_grassland 0.070347 0.048376 1.454 0.145903
pheno_endopen_forest 0.019023 0.090207 0.211 0.832977
pheno_endtree_grass 0.048338 0.060670 0.797 0.425597
sl_scaled -0.001803 0.008832 -0.204 0.838234
log_sl_scaled 0.154034 0.011179 13.779 < 2e-16
cos_ta_scaled 0.084257 0.007907 10.656 < 2e-16
sl_scaled:pheno_startbare_ground -0.009350 0.018311 -0.511 0.609606
sl_scaled:pheno_startdry_grassland 0.007096 0.021411 0.331 0.740318
sl_scaled:pheno_startopen_forest -0.090468 0.022801 -3.968 7.25e-05
sl_scaled:pheno_starttree_grass -0.040275 0.025024 -1.609 0.107510
log_sl_scaled:pheno_startbare_ground -0.074320 0.019119 -3.887 0.000101
log_sl_scaled:pheno_startdry_grassland -0.051027 0.025076 -2.035 0.041866
log_sl_scaled:pheno_startopen_forest -0.114017 0.023049 -4.947 7.55e-07
log_sl_scaled:pheno_starttree_grass -0.156630 0.024821 -6.310 2.78e-10
cos_ta_scaled:pheno_startbare_ground -0.104274 0.014004 -7.446 9.63e-14
cos_ta_scaled:pheno_startdry_grassland -0.013069 0.017470 -0.748 0.454429
cos_ta_scaled:pheno_startopen_forest -0.156955 0.017241 -9.104 < 2e-16
cos_ta_scaled:pheno_starttree_grass -0.141402 0.018059 -7.830 4.89e-15
(Intercept)
pheno_endbare_ground ***
pheno_enddry_grassland
pheno_endopen_forest
pheno_endtree_grass
sl_scaled
log_sl_scaled ***
cos_ta_scaled ***
sl_scaled:pheno_startbare_ground
sl_scaled:pheno_startdry_grassland
sl_scaled:pheno_startopen_forest ***
sl_scaled:pheno_starttree_grass
log_sl_scaled:pheno_startbare_ground ***
log_sl_scaled:pheno_startdry_grassland *
log_sl_scaled:pheno_startopen_forest ***
log_sl_scaled:pheno_starttree_grass ***
cos_ta_scaled:pheno_startbare_ground ***
cos_ta_scaled:pheno_startdry_grassland
cos_ta_scaled:pheno_startopen_forest ***
cos_ta_scaled:pheno_starttree_grass ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 182.413 sec elapsedCode
Code
[1] "Model location: buffalo_pheno_mixed_diag_mvmt_10rs_shrubgrassref_wet_2025-05-20.rds"Plot both models
Create manuscript figures
Code
coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
"Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass",
"sl", "log(sl)", "cos(ta)",
"sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
"log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
"cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
),
"Estimate" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"],
"SE" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Std. Error"]
) %>%
mutate(
# Calculate confidence intervals
LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE,
UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
# Compute p-values
pvalue = 2 * pnorm(-abs(Estimate) / SE)
)
# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
if (coef_df$pvalue[x] <= 0.001){
return("***")
} else if (coef_df$pvalue[x] <= 0.01) {
return("**")
} else if (coef_df$pvalue[x] <= 0.05) {
return("*")
}
})
# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")
# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))
# Specify the order in which the coefficients should be plotted
order <- c(
"Floodplain",
"DryGrassland",
"OpenForest",
"TreeGrass",
"sl",
"log(sl)",
"cos(ta)",
"sl:Floodplain",
"sl:DryGrassland",
"sl:OpenForest",
"sl:TreeGrass",
"log(sl):Floodplain",
"log(sl):DryGrassland",
"log(sl):OpenForest",
"log(sl):TreeGrass",
"cos(ta):Floodplain",
"cos(ta):DryGrassland",
"cos(ta):OpenForest",
"cos(ta):TreeGrass"
)
# Specify the order in which the coefficients should be plotted
order_no_mvmt <- c(
"Floodplain",
"DryGrassland",
"OpenForest",
"TreeGrass"
)
# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))
coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))Plot the coefficients
Code
min_x <- -0.5
max_x <- 0.8
# coefficients on the x-axis
dry_mvmt_plot <- ggplot(data = coef_df,
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 12.5, yend = 12.5,
colour = "gray80", lty = 1, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 15.5, yend = 15.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
# drop the filter argument when plotting all covariates
geom_errorbarh(data = coef_df2,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#FF8123", "#9B4200")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
# ggtitle("Dry season") +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#FF8123")))
dry_mvmt_plotWarning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 9 rows containing missing values or values outside the scale range
(`geom_text()`).Drop the filter argument when plotting all covariates
Code
dry_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")),
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
# drop the filter argument when plotting all covariates
geom_errorbarh(data = coef_df2_no_mvmt,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
# scale_y_discrete(limits = rev(order)) +
scale_y_discrete(limits = rev(order_no_mvmt)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#FF8123", "#9B4200")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#FF8123")))
dry_mvmt_plotWarning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).Code
coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
"Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass",
"sl", "log(sl)", "cos(ta)",
"sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
"log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
"cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
),
"Estimate" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"],
"SE" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Std. Error"]
) %>%
mutate(
# Calculate confidence intervals
LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE,
UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
# Compute p-values
pvalue = 2 * pnorm(-abs(Estimate) / SE)
)
# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
if (coef_df$pvalue[x] <= 0.001){
return("***")
} else if (coef_df$pvalue[x] <= 0.01) {
return("**")
} else if (coef_df$pvalue[x] <= 0.05) {
return("*")
}
})
# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")
# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))
# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))
coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
spread(key = Type, value = value) %>%
mutate(Level = paste0(Level, "%"))Plot the coefficients
Code
# coefficients on the x-axis
wet_mvmt_plot <- ggplot(data = coef_df,
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 12.5, yend = 12.5,
colour = "gray80", lty = 1, lwd = 0.3) +
annotate(geom = "segment",
x = min_x, xend = max_x,
y = 15.5, yend = 15.5,
colour = "gray80", lty = 1, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
geom_point(shape = 3, size = 2.5) +
geom_errorbarh(data = coef_df2,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
scale_y_discrete(limits = rev(order)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#238DFF", "#004691")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#1377E2")))
wet_mvmt_plotWarning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).Drop the filter argument when plotting all covariates
Code
wet_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")),
aes(y = Covariate, x = Estimate, col = factor(Preference))) +
geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
geom_point(shape = 3, size = 2.5) +
geom_errorbarh(data = coef_df2_no_mvmt,
aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)),
height = 0, alpha = 0.5) +
geom_text(aes(label = Significance, hjust = 0.5, vjust = 0),
show.legend = F) +
scale_y_discrete(limits = rev(order_no_mvmt)) +
theme_classic() +
scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.25)) +
coord_capped_cart(left = "both", bottom = "both") +
labs(x = expression(beta*"-Estimate")) +
scale_color_manual(name = "Preference",
values = c("#238DFF", "#004691")) +
scale_linewidth_manual(name = "Confidence Level",
values = c(2, 1, 0.3)) +
theme(legend.position = "bottom",
legend.margin = margin(0, 50, 0, -20),
legend.box.margin = margin(-5, -10, -5, -10),
legend.text = element_text(face = 3),
legend.title = element_text(face = 3)) +
guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
linewidth = guide_legend(title.position = "top", title.hjust = 0.5,
override.aes = list(colour = "#1377E2")))
wet_mvmt_plotWarning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).
Code
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).Updating movement parameters, with vegetation interactions
Code
coef_df_dry_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
"scaled_coefficient" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"])
# step length coefficients
sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate))
sl_dry_mvmt_mvparams <- sl_dry_mvmt_mvparams %>% mutate(
scale_factor = sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# log step length coefficients
log_sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate))
log_sl_dry_mvmt_mvparams <- log_sl_dry_mvmt_mvparams %>% mutate(
scale_factor = log_sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# turning angle coefficients
cos_ta_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate))
cos_ta_dry_mvmt_mvparams <- cos_ta_dry_mvmt_mvparams %>% mutate(
scale_factor = cos_ta_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
tentative_mean_sl <- tentative_shape * tentative_scale
# shrub_grass
shrub_grass_mvmt_dry <- data.frame(
"pheno" = "ShrubGrass",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1],
"updated_scale" = 1/((1/tentative_scale) - sl_dry_mvmt_mvparams$nat_coef[1]),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1]
)
shrub_grass_mean_sl_dry <- shrub_grass_mvmt_dry$updated_shape * shrub_grass_mvmt_dry$updated_scale
# floodplain
floodplain_mvmt_dry <- data.frame(
"pheno" = "Floodplain",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[2],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[2])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[2]
)
floodplain_mean_sl_dry <- floodplain_mvmt_dry$updated_shape * floodplain_mvmt_dry$updated_scale
# dry_grassland
dry_grassland_mvmt_dry <- data.frame(
"pheno" = "DryGrassland",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[3],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[3])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[3]
)
dry_grassland_mean_sl_dry <- dry_grassland_mvmt_dry$updated_shape * dry_grassland_mvmt_dry$updated_scale
# open_forest
open_forest_mvmt_dry <- data.frame(
"pheno" = "OpenForest",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[4],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[4])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[4]
)
open_forest_mean_sl_dry <- open_forest_mvmt_dry$updated_shape * open_forest_mvmt_dry$updated_scale
# tree_grass
tree_grass_mvmt_dry <- data.frame(
"pheno" = "TreeGrass",
"updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[5],
"updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[5])),
"updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[5]
)
tree_grass_mean_sl_dry <- tree_grass_mvmt_dry$updated_shape * tree_grass_mvmt_dry$updated_scaleCreating Gamma distributions
Code
sl_dry_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))
# Tentative distribution
sl_dry_mvmt$Tentative <- dgamma(x = sl_dry_mvmt$x,
shape = tentative_shape,
scale = tentative_scale)
# Shrub grass
sl_dry_mvmt$ShrubGrass <- dgamma(x = sl_dry_mvmt$x,
shape = shrub_grass_mvmt_dry$updated_shape,
scale = shrub_grass_mvmt_dry$updated_scale)
# Bare ground
sl_dry_mvmt$Floodplain <- dgamma(x = sl_dry_mvmt$x,
shape = floodplain_mvmt_dry$updated_shape,
scale = floodplain_mvmt_dry$updated_scale)
# Dry grassland
sl_dry_mvmt$DryGrassland <- dgamma(x = sl_dry_mvmt$x,
shape = dry_grassland_mvmt_dry$updated_shape,
scale = dry_grassland_mvmt_dry$updated_scale)
# Open forest
sl_dry_mvmt$OpenForest <- dgamma(x = sl_dry_mvmt$x,
shape = open_forest_mvmt_dry$updated_shape,
scale = open_forest_mvmt_dry$updated_scale)
# Tree grass
sl_dry_mvmt$TreeGrass <- dgamma(x = sl_dry_mvmt$x,
shape = tree_grass_mvmt_dry$updated_shape,
scale = tree_grass_mvmt_dry$updated_scale)
# Pivot from wide to long data
sl_dry_mvmt <- sl_dry_mvmt %>%
pivot_longer(cols = -x)
# Specify the order in which the coefficients should be plotted
order <- c(
"ShrubGrass",
"Floodplain",
"DryGrassland",
"OpenForest",
"TreeGrass",
"Tentative"
)
sl_dry_mvmt$name <- factor(sl_dry_mvmt$name, levels = order)
mean_values <- c("ShrubGrass" = round(shrub_grass_mean_sl_dry, 1),
"Floodplain" = round(floodplain_mean_sl_dry, 1),
"DryGrassland" = round(dry_grassland_mean_sl_dry, 1),
"OpenForest" = round(open_forest_mean_sl_dry, 1),
"TreeGrass" = round(tree_grass_mean_sl_dry, 1),
"Tentative" = round(tentative_mean_sl, 1))
# Construct dynamic labels
labels <- paste(names(mean_values), ": Mean SL = ", mean_values, sep = "")
# Plot
ggplot(sl_dry_mvmt, aes(x = x, y = value, color = name, linetype = name)) +
geom_line(alpha = 0.75, linewidth = 1) +
scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black"),
labels = labels) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted"),
labels = labels) +
theme_bw() +
theme(legend.position = c(0.65, 0.65))Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
Creating von Mises distributions
Code
ta_dry_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))
# Tentative distribution
ta_dry_mvmt$Tentative <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = tentative_kappa)
# Shrub grass
ta_dry_mvmt$ShrubGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = shrub_grass_mvmt_dry$updated_kappa)
# Bare ground
ta_dry_mvmt$Floodplain <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = floodplain_mvmt_dry$updated_kappa)
# Dry grassland
ta_dry_mvmt$DryGrassland <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = dry_grassland_mvmt_dry$updated_kappa)
# Open forest
ta_dry_mvmt$OpenForest <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = open_forest_mvmt_dry$updated_kappa)
# Tree grass
ta_dry_mvmt$TreeGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
mu = 0,
kappa = tree_grass_mvmt_dry$updated_kappa)
# Pivot from wide to long data
ta_dry_mvmt <- ta_dry_mvmt %>%
pivot_longer(cols = -x)
ta_dry_mvmt$name <- factor(ta_dry_mvmt$name, levels = order)
# Plot
ggplot(ta_dry_mvmt,
aes(x = x, y = value, color = factor(name), linetype = factor(name))
) +
geom_line(alpha = 0.75, linewidth = 1) +
xlab("Turning angle (radians)") +
scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black")) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted")) +
# scale_colour_viridis_d("Habitat", option = "D") +
# ggtitle("Dry season") +
theme_bw() +
theme(legend.position = "bottom")
Code
coef_df_wet_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
"scaled_coefficient" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"])
# step length coefficients
sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate))
sl_wet_mvmt_mvparams <- sl_wet_mvmt_mvparams %>% mutate(
scale_factor = sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# log step length coefficients
log_sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate))
log_sl_wet_mvmt_mvparams <- log_sl_wet_mvmt_mvparams %>% mutate(
scale_factor = log_sl_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# turning angle coefficients
cos_ta_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate))
cos_ta_wet_mvmt_mvparams <- cos_ta_wet_mvmt_mvparams %>% mutate(
scale_factor = cos_ta_scale_sd,
nat_coef = scaled_coefficient / scale_factor,
pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
# shrub_grass
shrub_grass_mvmt_wet <- data.frame(
"pheno" = "ShrubGrass",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1],
"updated_scale" = 1/((1/tentative_scale) - sl_wet_mvmt_mvparams$nat_coef[1]),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1]
)
shrub_grass_mean_sl_wet <- shrub_grass_mvmt_wet$updated_shape * shrub_grass_mvmt_wet$updated_scale
# floodplain
floodplain_mvmt_wet <- data.frame(
"pheno" = "Floodplain",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[2],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[2])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[2]
)
floodplain_mean_sl_wet <- floodplain_mvmt_wet$updated_shape * floodplain_mvmt_wet$updated_scale
# dry_grassland
dry_grassland_mvmt_wet <- data.frame(
"pheno" = "wetGrassland",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[3],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[3])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[3]
)
dry_grassland_mean_sl_wet <- dry_grassland_mvmt_wet$updated_shape * dry_grassland_mvmt_wet$updated_scale
# open_forest
open_forest_mvmt_wet <- data.frame(
"pheno" = "OpenForest",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[4],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[4])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[4]
)
open_forest_mean_sl_wet <- open_forest_mvmt_wet$updated_shape * open_forest_mvmt_wet$updated_scale
# tree_grass
tree_grass_mvmt_wet <- data.frame(
"pheno" = "TreeGrass",
"updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[5],
"updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[5])),
"updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[5]
)
tree_grass_mean_sl_wet <- tree_grass_mvmt_wet$updated_shape * tree_grass_mvmt_wet$updated_scaleCreating Gamma distributions
Code
sl_wet_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))
# Tentative distribution
sl_wet_mvmt$Tentative <- dgamma(x = sl_wet_mvmt$x,
shape = tentative_shape,
scale = tentative_scale)
# Shrub grass
sl_wet_mvmt$ShrubGrass <- dgamma(x = sl_wet_mvmt$x,
shape = shrub_grass_mvmt_wet$updated_shape,
scale = shrub_grass_mvmt_wet$updated_scale)
# Bare ground
sl_wet_mvmt$Floodplain <- dgamma(x = sl_wet_mvmt$x,
shape = floodplain_mvmt_wet$updated_shape,
scale = floodplain_mvmt_wet$updated_scale)
# Dry grassland
sl_wet_mvmt$DryGrassland <- dgamma(x = sl_wet_mvmt$x,
shape = dry_grassland_mvmt_wet$updated_shape,
scale = dry_grassland_mvmt_wet$updated_scale)
# Open forest
sl_wet_mvmt$OpenForest <- dgamma(x = sl_wet_mvmt$x,
shape = open_forest_mvmt_wet$updated_shape,
scale = open_forest_mvmt_wet$updated_scale)
# Tree grass
sl_wet_mvmt$TreeGrass <- dgamma(x = sl_wet_mvmt$x,
shape = tree_grass_mvmt_wet$updated_shape,
scale = tree_grass_mvmt_wet$updated_scale)
# Pivot from wide to long data
sl_wet_mvmt <- sl_wet_mvmt %>%
pivot_longer(cols = -x)
sl_wet_mvmt$name <- factor(sl_wet_mvmt$name, levels = order)
mean_values_wet <- c("ShrubGrass" = round(shrub_grass_mean_sl_wet, 1),
"Floodplain" = round(floodplain_mean_sl_wet, 1),
"DryGrassland" = round(dry_grassland_mean_sl_wet, 1),
"OpenForest" = round(open_forest_mean_sl_wet, 1),
"TreeGrass" = round(tree_grass_mean_sl_wet, 1),
"Tentative" = round(tentative_mean_sl, 1))
# Construct dynamic labels
labels <- paste(names(mean_values_wet), ": Mean SL = ", mean_values_wet, sep = "")
# Plot
ggplot(sl_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
geom_line(alpha = 0.75, linewidth = 1) +
scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black"),
labels = labels) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted"),
labels = labels) +
# scale_colour_viridis_d("Habitat", option = "D") +
# ggtitle("Wet season") +
theme_bw() +
theme(legend.position = c(0.65, 0.65))
Creating von Mises distributions
Code
ta_wet_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))
# Tentative distribution
ta_wet_mvmt$Tentative <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = tentative_kappa)
# Shrub grass
ta_wet_mvmt$ShrubGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = shrub_grass_mvmt_wet$updated_kappa)
# Bare ground
ta_wet_mvmt$Floodplain <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = floodplain_mvmt_wet$updated_kappa)
# Dry grassland
ta_wet_mvmt$DryGrassland <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = dry_grassland_mvmt_wet$updated_kappa)
# Open forest
ta_wet_mvmt$OpenForest <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = open_forest_mvmt_wet$updated_kappa)
# Tree grass
ta_wet_mvmt$TreeGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
mu = 0,
kappa = tree_grass_mvmt_wet$updated_kappa)
# Pivot from wide to long data
ta_wet_mvmt <- ta_wet_mvmt %>%
pivot_longer(cols = -x)
ta_wet_mvmt$name <- factor(ta_wet_mvmt$name, levels = order)
# Plot
ggplot(ta_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
geom_line(alpha = 0.75, linewidth = 1) +
xlab("Turning angle (radians)") +
scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
scale_colour_manual("Habitat",
values = c("ShrubGrass" = "#999900",
"Floodplain" = "#3399FF",
"DryGrassland" = "#CC9900",
"OpenForest" = "#006600",
"TreeGrass" = "#00CC66",
"Tentative" = "black")) +
scale_linetype_manual("Habitat",
values = c("ShrubGrass" = "solid",
"Floodplain" = "solid",
"DryGrassland" = "solid",
"OpenForest" = "solid",
"TreeGrass" = "solid",
"Tentative" = "dotted")) +
# scale_colour_viridis_d("Habitat", option = "D") +
# ggtitle("Dry season") +
theme_bw() +
theme(legend.position = "bottom")
Session Information
R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.4.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Australia/Brisbane
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] amt_0.2.2.0 circular_0.5-1 lemon_0.5.0
[4] beepr_2.0 tictoc_1.2.1 glmmTMB_1.1.11
[7] sjPlot_2.8.17 TwoStepCLogit_1.2.6 lubridate_1.9.4
[10] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[13] purrr_1.0.4 readr_2.1.5 tidyr_1.3.1
[16] tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] sjlabelled_1.2.0 tidyselect_1.2.1 farver_2.1.2
[4] fastmap_1.2.0 bayestestR_0.15.3 sjstats_0.19.0
[7] digest_0.6.37 timechange_0.3.0 lifecycle_1.0.4
[10] sf_1.0-20 survival_3.8-3 magrittr_2.0.3
[13] compiler_4.5.0 rlang_1.1.6 tools_4.5.0
[16] yaml_2.3.10 knitr_1.50 labeling_0.4.3
[19] bit_4.6.0 classInt_0.4-11 plyr_1.8.9
[22] RColorBrewer_1.1-3 KernSmooth_2.23-26 withr_3.0.2
[25] numDeriv_2016.8-1.1 grid_4.5.0 datawizard_1.1.0
[28] e1071_1.7-16 scales_1.4.0 MASS_7.3-65
[31] insight_1.2.0 cli_3.6.5 mvtnorm_1.3-3
[34] crayon_1.5.3 rmarkdown_2.29 ragg_1.4.0
[37] reformulas_0.4.1 generics_0.1.3 rstudioapi_0.17.1
[40] performance_0.13.0 tzdb_0.5.0 parameters_0.25.0
[43] minqa_1.2.8 DBI_1.2.3 proxy_0.4-27
[46] audio_0.1-11 splines_4.5.0 parallel_4.5.0
[49] effectsize_1.0.0 vctrs_0.6.5 boot_1.3-31
[52] Matrix_1.7-3 jsonlite_2.0.0 hms_1.1.3
[55] bit64_4.6.0-1 systemfonts_1.2.3 units_0.8-7
[58] glue_1.8.0 nloptr_2.2.1 stringi_1.8.7
[61] gtable_0.3.6 ggeffects_2.2.1 lme4_1.1-37
[64] pillar_1.10.2 htmltools_0.5.8.1 R6_2.6.1
[67] TMB_1.9.17 textshaping_1.0.1 Rdpack_2.6.4
[70] vroom_1.6.5 evaluate_1.0.3 lattice_0.22-6
[73] rbibutils_2.3 class_7.3-23 Rcpp_1.0.14
[76] gridExtra_2.3 nlme_3.1-168 mgcv_1.9-1
[79] xfun_0.52 sjmisc_2.8.10 pkgconfig_2.0.3